Near mean-field behavior in the generalized Burridge-Knopoff 
earthquake model with variable range stress transfer 



J3 



Junchao Xia_| and Harvey Gould 
Department of Physics, Clark University, Worcester, MA 01610 

W. Klein 

Department of Physics and Center for Computational Science, 
Boston University, Boston, MA 02215 

J. B. Rundle 

Department of Physics and Center for Computational Science and Engineering, 
University of California, Davis, CA 95616 

Simple models of earthquake faults are important for understanding the mecha- 
nisms for their observed behavior in nature, such as Gutenberg-Richter scaling. Be- 
cause of the importance of long-range interactions in an elastic medium, we general- 
ize the Burridge-Knopoff slider-block model to include variable range stress transfer. 
We find that the Burridge-Knopoff model with long-range stress transfer exhibits 
qualitatively different behavior than the corresponding long-range cellular automata 
models and the usual Burridge-Knopoff model with nearest-neighbor stress transfer, 
depending on how quickly the friction force weakens with increasing velocity. Exten- 
sive simulations of quasiperiodic characteristic events, mode-switching phenomena, 
ergodicity, and waiting-time distributions are also discussed. Our results are consis- 
tent with the existence of a mean-field critical point and have important implications 
for our understanding of earthquakes and other driven dissipative systems. 
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I. INTRODUCTION 



Earthquake faults are important examples of driven dissipative systems Models of 
fault systems are important for understanding Gutenberg- Richter (power law) behavior, the 
occurrence of characteristic events, and the relation between small and large earthquakes l|, 
Understanding driven dissipative systems is important, for example, in the 
I nenra. networks g depinning transition, in charge density waves 
magnetized domains in ferromagnets M , domain rearrangements 
in flowing foams 10], and granular materials under shear stress [111 ]. 

A relatively simple dynamical model that contains much of the essential physics of earth- 



context of avalanches 
and superconductors 



quake faults is the spring-block model proposed by Burridge and Knopoff [12]. This model 



consists of blocks connected by linear springs to their nearest neighbors with spring constant 
k c . The blocks are also connected to a loader plate by linear springs with spring constant 
ki, and rest on a surface with a nonlinear velocity- weakening stick-slip friction force which 
depends on a parameter a that controls how quickly the friction force decreases as the ve- 
locity is increased. The model was studied numerically in one dimension in Ref. [3] and 



more recently in Refs. 



0, Q, H, 
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17. 
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An earthquake event is defined as a cluster of blocks that move (slip) due to the initial 
slip of a single block. In addition to the amount of energy released in an earthquake event, 
a quantity of interest is the moment M, which is defined as . Axj, where Axj is the net 
displacement of block j during an event and the sum is over all the blocks in the event. 
Carlson and Langer simulated the one-dimensional Burridge-Knopoff model for N = 100 
and iV = 1000 blocks. The main result of their simulations 1 1 



14 



ipo : 

y, 



16 



113, llSj is that for 

a > 2 the moment probability distribution P(M) scales as M~ x for small localized events 
with an exponent x ~ 2 j^]. There also is a peak in P(M) for large events indicating a 
significant presence of characteristic (non-power law) events. 

Because simulations of the Burridge-Knopoff model require solving Newton's equations 
of motion and are time consuming, several cellular automata (CA) models have been pro- 
posed that neglect the inertia of the blocks and simplify the effect of the friction force by 
assuming that the motion is overdamped. These cellular automata include those due to 



Rundle, Jackson, and Brown 



23] and Olami, Feder, and Christensen 24J. In these models 



P(s), the distribution of the number of blocks in an event, does not exhibit power law seal- 
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mg 



25] for nearest-neighbor stress transfer if periodic boundary conditions are used 26) . A 



generalization of these CA models to include more realistic long-range stress transfer 27] 
yields considerable differences with the nearest-neighbor CA models 28j and with the orig- 
inal Burridge-Knopoff model. In particular, for long-range stress transfer P(s) exhibits 
Gutenberg-Richter scaling consistent with the system being near a mean-field (spinodal) 



critical point 



29 



30, 



a Langevin equation 



31 

29|, 



31 



. In addition, the long-range CA models can be described by 



331 ] and small and medium size events can be interpreted as 



fluctuations about a free energy minimum [29j, |3Cj, |34|]. Large events drive the system out of 
equilibrium from which the system decays back to an equilibrium state {29]]. 

The CA and Burridge-Knopoff models lack several elements that would make them more 
realistic. In particular, the long-range CA models do not include inertia and more realistic 
friction laws, and the Burridge-Knopoff model does not include long-range stress transfer. 
Both types of models do not include elastic (seismic) wave radiation (phonons) because there 
is no medium in which seismic waves can propagate. However, the lack of seismic waves is 
a reasonable approximation, because seismic waves carry little energy in real faults 35]. 

In this paper we discuss our extensive simulations of a generalized Burridge-Knopoff 



model with long-range stress transfer between the blocks [361 ] ■ For various values of the 
dynamic friction parameter a and the range of stress transfer R, we observed phenomena 
similar to real fault network systems, including Gutenberg-Richter scaling, quasiperiodic 
characteristic events, and mode-switching. Our primary results are that the behavior of 
the long-range Burridge-Knopoff model differs significantly from the short-range Burridge- 
Knopoff model, the behavior of the long-range Burridge-Knopoff and CA models is similar 
only for small a, and the nature of the friction force is important and strongly affects the 
behavior of the Burridge-Knopoff model. In particular, we find numerical evidence for two 

d spinodal critical point similar to that found in the 



types of scaling behavior: a mean-fie 
lone-range CA models 

found in Refs. Q, [l5, 19] for a ^ 2 and all values of R studied in the range 1 < R < 500. 



30 



31 



32] for R 3> 1 and a ^ 1 36] and the scaling behavior 
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II. BURRIDGE-KNOPOFF MODEL 



The orig inal Burridge-Knopoff model in one dimension is governed by the equation of 



motion 112 



m- 



dt 2 



k c (xj+i — 2xj + Xj-\) — kiXj — F{y + Xj), 



(1) 



where Xj is the displacement of block j from its equilibrium position, v is the speed of the 
substrate, which moves to the left, F(x) = Fq4>(x/v) is a velocity-dependent friction force, 
v is a characteristic velocity, and m is the mass of a block. The loader plate is fixed. 

- - k L /m, and Uj = (k L /F )xj, 



As in Ref. [ijj] we introduce the scaled variables r 
and rewrite Eq. (pQ) in dimensionless form as 



Uj = £ (uj + i — 2uj + — Uj — (f)(2av + 2aUj) : 



(2) 



with 2a = u p F /ki,v, 



k c /k L , and v = vkL/(u} p F ); a dot denotes differentiation with 



respect to r. The form of the friction force is plotted in Fig. [T]and is given by 



-oo, 1], 
1 - a 



1 + 



y 



y = 
y>0. 



(3) 



1 - a 



Note that <p(y) decays monotonically to zero from 4>(0 + ) = 1 — a and prohibits slip in the 
same direction as the motion of the substrate [371 ] . 

The four dimensionless parameters, £, a, u, and a govern the behavior of the system. 
The parameter a appears in the argument of (f) in Eq. (j2j) and determines how quickly the 
dynamic friction force decreases with increasing velocity; a = means that the dynamic 
friction force is equal to the constant 1 — a. Larger a means that the friction force decreases 
more rapidly with velocity and the motion is less damped; a — *■ oo implies that the friction 
force drops to zero immediately for positive velocities. 

We generalize the Burridge-Knopoff model by assuming that a block is connected to R 
neighbors (in each direction) with the rescaled spring constant k c j R 38] ; R = 1 corresponds 



to the usua' 



Burridge-Knopoff model. We used the second- and fourth-order Runge-Kutta 



algorithms 39|, |40j with the time step At = 0.001 to solve Eq. (j2J) generalized to arbitrary 



R. Both algorithms and other fourth-order algorithms 



4l| give similar results. 
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FIG. 1: The form of the velocity- weakening friction force <j)(y). The friction force decays monoton- 
ically to zero from the initial value (j)(0 + ) = 1 — a with the initial slope —1. 

III. IMPLEMENTATION OF THE BURRIDGE-KNOPOFF MODEL 

Because the velocity of a block is a continuous variable, we need to introduce a criterion 
for when a "stuck" block begins to move and when a moving block becomes stuck so that we 
can define the beginning and end of an event. We define a block to be stuck if its velocity 
is less than a parameter vq. In addition, the stress on the block, defined to be the force due 
to all the springs coupled to it including the loader plate spring, must be smaller than the 
maximum static friction force Fq (taken to be unity in dimensionless units). If a block is 
stuck, we choose the value of the static friction force to be such that it cancels the stress. 
At the next time step, a stuck block will remain stuck if the stress on it is still smaller than 
Fq. A moving block will become stuck at the next time step if its speed is less than vq and 
decreasing and if the stress on it is less than Fq. In our simulations we take vq = 10~ 5 , which 
yields reasonable results. 

An earthquake event begins with the slip of a block and ends when all blocks become 
stuck. A block is said to "fail" when it begins to move after being stuck. A moving block 
can become stuck and then move (slip) again during an event. 

We initially set itj = for all j and assign small random displacements to all the blocks; 
hence all blocks are initially stuck. We compute the force on all the blocks and update 
itj and Uj for all j using the generalization of Eq. (j2j) for arbitrary R. We continue these 
updates until all blocks become stuck again. We then move the substrate (the loader plate 
is fixed) until the stress on one block exceeds Fq. This stress loading mechanism is known 
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FIG. 2: (Color online) Log-log plot of P(M) versus M for different values of a for R = 1. Scaling 
behavior can be found for small events for a c 1. The exponent x defined by P(M) ~ M _x 
decreases from x ~ 2 as a is decreased from 2.5 to 1. For a < 1, no scaling behavior is observed. 
A nonuniform bin size was used here and in Figs. [3] and [H For clarity each distribution is shifted 
vertically by 2 units. 



as the zero velocity limit 15j, |32|, |42j and is equivalent to setting v = in Eq. The 



zero velocity limit is realistic because the dynamics of earthquake faults can be divided into 
continuous loading on a long time scale and relaxation with release of stress and energy 
on a much shorter time scale. The zero velocity limit ensures that there is only one event 
per substrate update and saves considerable simulation time. The relaxation process occurs 
when the motion or failure of the initiator induces other blocks to slip. If there are "moving" 
blocks, the event is still alive; otherwise, we reload the system to induce a new event. 

The results in this paper are for a = 0.01, £ = 10, Vq = 10~ 5 , iV = 5000, 10 6 events, 
and various values of R and a. Most previous work has been for R = 1 with N = 100 [ijj], 



N = 1000 N = 800 



20 



are used as in previous work 14 



2l| . and the same values of o and £. Open boundary conditions 



15 



20. 
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IV. THE SIZE DISTRIBUTION OF EARTHQUAKE EVENTS 

In Refs. [14! [15I the moment distribution P(M) was found to exhibit small localized 



events and larger delocalized events (for M > 21/ a). The localized events show power-law 
behavior with a slope x ~ 2; the delocalized or characteristic earthquakes correspond to a 
pronounced peak in P(M). Carlson and Langer [15I considered £ = 6, 8, 10, and 12 
and a = 2.5, 3, 4, and 5, and found the same general behavior for P(M). They also found 
that the effective value of the exponent x decreases as a — > 1 and the power law behavior 
becomes less well defined. 

Figure [2] shows our results for the moment distribution P(M) for R = 1 and values of a 
in the range 0.5 < a < 2.5. We find that for a ?t 1, P(M) exhibits power law behavior with 
slope x « 2 for small localized events in the range 10~ 4 < M ^ 10° and a non-power law 
distribution of characteristic events. For a < 1, no power law behavior is found [43j. These 



14,05m 



results are consistent with the results in Refs. 

Figure [3] shows P(M) for R = 500 and the same values of a as in Fig. [2j For a > 1 we 
see that the increased interaction range R does not change the value of the exponent x, but 
the scaling range becomes smaller, and the distribution of characteristic events (M ^ 1) is 
concentrated in a narrower range close to the system size. We conclude that although P(M) 
exhibits power law behavior for a = 2.5, the slope x ~ 2 differs from the mean-field value 



of x = 3/2 found in the long-range CA models [3l|, |32 |. 

In contrast, for R > 100 and a ^ 1 the scaling exponent approaches x m 1.5 as a — » 0. 
This value of x is close to the mean-field value of 3/2 found in the long-range CA models. 
Increasing R increases the range of the power law behavior and decreases the number of 
characteristic events. For a = the mean-field behavior of P(M) becomes even better 
defined. We conclude that the scaling exponent of P(M) converges to the mean-field value 
of 3/2 for sufficiently small values of a and sufficiently large values of R. 

To compare our results more directly with the cellular automata models, we compute 
P{s), the distribution of the number of blocks in an event (including multiple failures). As 



shown in Fig. 4(a) , there is no power law behavior for R = 1 and all values of a studied 
in contrast to the power law behavior of P(M). This result is similar to that observed in 
Ref. 15j. Our results for P(s) with R = 500 are shown in Fig. 4(b) for the same values of a 



as in Fig. 4(a) We see that for R = 500, P(M) and P(s) display similar power law behavior 
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with the slope x ~ 2 for a ~ 2.5 and x m 1.5 for a < 1 (see Tabled]). As expected, the power 
law behavior of P(s) and P(M) does not hold for large s. These large scale characteristic 
events are also observed in the cellular automata models and become less probable as R is 



increased 



44 



45 



46]. 



The scaling behavior for a = 3, 4, 5 and 10 is similar to that for a = 2.5. We conclude that 
the generalized Burridge-Knopoff model exhibits Gutenberg-Richter scaling with a mean- 
field exponent of 1.5 for small a and large R. Different scaling behavior is found for a > 2 
and all values of R studied. 



V. THE MEAN SLIP, STRESS DROP, AND THE NUMBER OF FAILURES 



In the CA models [23, 



24 



25 



31 



, [32J the stress on a block after it fails decreases to its 



residual stress, which is chosen to be either the same for all blocks, or if noise is introduced, 



Oh 

O 

(50 

o 



-10 



-15 



-20 



a+= 2.5 

+ 

0k=2.0 + +. 

x + ++ ++ 
<M= 1.5 x x x ++ + +J 

X 

Cb= 1.0 * 

□ 

Of = 0.5 D n n 
ot>= 0.0 ' 

o „ 




-4-3-2-101234 

logjo M 

FIG. 3: (Color online) Log-log plot of P(M) versus M for R = 500 and the same values of a as in 
Fig. [21 For a <j 1 the power law exponent converges to the mean- field value of x = 1.5 as a — ► 0. 
For a ^ 1, the scaling exponent becomes close to 2 as a is increased. The scaling range becomes 
smaller as a increases. For clarity each distribution is shifted vertically by 2 units. (The apparent 
linear behavior for M ^ 1 is an artifact due to limited statistics and the use of a nonlinear bin 
width.) 
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FIG. 4: (Color online) The distribution of events with s blocks, P(s). (a) Log- log plot of P(s) 
versus s for R = 1 and various values of a. In contrast to the behavior of P(M), P(s) does not 
exhibit power law behavior for R = 1 for all values of a studied, (b) For R = 500, P(s) ~ s~ x 
with x ~ 2 for a ~ 2.5 and x ~ 1.5 for a < 0.5. Note that s can exceed iV because a block can 
slip and become stuck and then slip again during an event. For clarity each distribution is shifted 
vertically by 2 units. 



is the same on the average. In the Rundle- Jackson-Brown model 23|, each failed block is 
displaced by an amount corresponding to the decrease of its stress. For large R the stress 
on a block at failure approaches the failure threshold and therefore all blocks in an event 
are displaced the same amount. Hence, the moment M of an event and s, the total num- 



ber of 
model 



ai 



ed blocks in an event, are proportional for the long-range Rundle- Jackson-Brown 



3l| . and the scaling behavior of P(M) and P(s) becomes identical. In addition, in the 



mean-field limit, a site fails (slips) only once during an event (3l|, |32j. These conditions are 
assumed by the coarse-graining theory of the long-range Rundle- Jackson-Brown model 31] 
and have been verified by computer simulations [32]. We check here if these assumptions 
hold for the Burridge-Knopoff model for sufficiently large R. 

Figure [5] shows <Aw>, the mean displacement or slip of a block during an event, as a 
function of s. Each block is counted only once even if it slips multiple times. The behavior 
of <Aw> as a function of s is similar for a = 2.5 and a = 0.5, and the range of s over which 
<Au> is independent of s increases with the interaction range R for all values of a studied. 
The implication of this independence is that each block slips the same amount. Note that 
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no scaling 
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random 
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no scaling 


yes 


random 
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no scaling 


no scaling 


yes 


random 
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2.02 


2.01 


no 


quasi-periodic 
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1.82 


no 


quasi-periodic 


1 


500 


1.63 


1.62 


no 


quasi-periodic 


0.9 


500 


1.62 


1.60 


no 


mode-switching 



TABLE I: Summary of the behavior of the Burridge-Knopoff model for several values of R and a. 
For a > 1, P(M) exhibits power law behavior for all R studied with the exponent x m ~ 2. Power 
law behavior of P(s) is found only for R > 100; the corresponding exponent is denoted as x s . For 
a < 1, P(M) and P{s) exhibit mean-field scaling with slope x ~ 1.5 for R > 100. Ergodicity is 
discussed in Sec. Wl\ and the time series of the stress is discussed in Sec. IVII1 

for R = 1 there is no range of s over which <Aw> is independent of s, even though P(M) 
exhibits power law behavior for small M and a = 2.5. Also note that the mean slip of a 
block in a characteristic event increases with a for fixed R. 

Figure displays the moment M as a function of s. As expected, the range of s for which 
Mas increases with R because the range of s for which the mean displacement is constant 
increases with R. 

The mean decrease of the stress on a block after an event is given by 

<Af> = - s J2^h (4) 

i 

where the sum is over all blocks in the event, and A/j is the difference of the stress on the 
ith block before and after an event. In Fig. [7]we see that the range of constant <Af> scales 
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FIG. 5: (Color online) Log-log plot of the mean displacement (slip) <Au> of a block as a function 
of the number of blocks s in an event for (a) a = 2.5 and (b) a = 0.5. Each block is counted only 
once even if it fails multiple times. The mean slip of a block becomes independent of s over a wide 
range of s for R 3> 1, and the range scales with R. For R S> 1, a characteristic event for a = 2.5 
has a much larger mean slip than that for a = 0.5. 
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FIG. 6: (Color online) Log-log plot of the moment Mas a function of s, the number of blocks 
in an event for (a) a = 2.5 and (b) a = 0.5. The range of s for which M a s increases as R is 
increased. A log-log plot is used to show the linear region more clearly. 

with R; <Af> ?a a = 0.01 for power law (small s) events, independent of the value of a. 

The independence of the mean displacement and mean stress drop on the size of an earth- 
quake has been observed in real earthquakes and has been interpreted as evidence for their 
self-similarity . This independence on the size of an earthquake is one of the assumptions 
of the static crack model of earthquakes and suggests that the Burridge-Knopoff model 
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with long-range stress transfer can capture some of the aspects of real earthquakes. 

Figure [8] shows the mean number of times a block fails during an event as a function 
of s. We see that the range of s over which a block fails only once scales with R; that is, 
multiple failures occur only for nonscaling events that are larger than 2R, the total number 
of neighbors of a block. This behavior is independent of a. We conclude that in the limit 
R — > oo, there are no blocks with multiple failures, consistent with the assumption made 
in the coarse-graining description of the Rundle- Jackson-Brown model in the mean-field 
limit 
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VI. ERGODICITY AND THE STRESS METRIC 



We characterize the nature of ergodicity in the Burridge-Knopoff model by determining 
the metric of the stress [47] . We take fi(t) to be a quantity associated with block i 

and define 



fi{t) = \fm')dt' 
i N - 

i=l 



(5) 
(6) 
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FIG. 7: (Color online) The mean stress drop per block, <A/>, as a function of s for (a) a = 2.5 
and (b) a = 0.5. The decrease is approximately equal to a = 0.01 for small s. The range of s over 
which the decrease is independent of s is proportional to the interaction range R. Note that the 
stress drop in a characteristic event for a = 2.5 is a much larger than for a = 0.5. 
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FIG. 8: (Color online) The mean number of times a block fails during an event as a function of s, 
the size of an event, for (a) a = 2.5 and (b) a = 0.5. We see that the range of s over which blocks 
fail only once is proportional to the interaction range R. Note that the blocks in a characteristic 
event fail many more times for a = 2.5 than for a = 0.5. The curves extrapolate to 1 as log s — > 0. 



and the metric 



the stress metric 



1 N 

n /(*) = ^E [/<(*) -</(*)>]' 

0- ; 



(7) 



If the system is ergodic, fi/(t) oc l/t Because the metric studied in the CA models is 



32], we choose fj to be the stress on block j just after an event. 
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FIG. 9: (Color online) The inverse stress metric versus the number of substrate updates (loading 
times) for R = 1. The system is ergodic for all the values of a studied, and the slope is an increasing 
function of a (see Table U for more values of a). 
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FIG. 10: (Color online) The normalized inverse stress metric versus the number of substrate updates 
(loading times) for R = 100 and 500 and (a) a = 2.5, (b) a = 0.5, and (c) a = 0. Note the different 
vertical scales. The system is nonergodic for larger values of R for a = 2.5. For a = 0.5 and a = 0, 
the slope of Q(0)/fl(t) becomes larger as R is increased. The system exhibits punctuated ergodicity 



for a = 0.5 and R = 100, similar to the behavior of the long-range CA models 
for more values of a and R. 



32]. See Ref. 



46] 



In Fig.[9]we show the normalized inverse stress metric Qf(0)/Qf(t) for R = 1 and different 
values of a. We see that Qf(0)/Qf(t) increases linearly with t. The mixing time r can be 
defined by the relation Qf(0)/flf(t) = t/r. We see that the mixing time r decreases with 
increasing a {^J. We conclude that the nearest-neighbor (R = 1) Burridge-Knopoff model 
is ergodic for all values of a studied. 

The system exhibits qualitatively different behavior for larger values of R. For a = 2.5 
the system is nonergodic for R > 100 during our observation time. In contrast, for a ^ 1 the 
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system remains ergodic as R is increased and the mixing time r decreases with increasing 



a as for R = 1 (see Figs. 10(b) and 10(c)). Note that for a = 0.5 and R = 100, the system 



displays punctuated ergodicity during our observation time, similar to the behavior of the 
long-range CA models [321 ] . 

Punctuated ergodicity has been observed in the Southern California fault system for some 



coarse graining conditions 



481 ] . Ergodicity can also be recovered for a simple, far- from 



equilibrium system with underlying chaotic dynamics serving as a temperature bath 
Our results suggest the possibility of both ergodic faults and nonergodic faults. 



49] 



VII. THE TIME SERIES OF THE STRESS 



To understand the behavior of the stress metric, we plot the times series of the stress / 
per block just after an event (see Fig. [TTl) . For a = 2.5 and R — 1, the times series fluctuates 
between to 0.4. This behavior is consistent with the observed ergodicity of the system (see 



Fig. [9]). As R is increased (Figjll(b) ), the time series becomes quasiperiodic and the period 



becomes longer for larger a 46j. For R = 500, f(t) is quasiperiodic with a mean of ~ 0.4 
and a range from « —0.2 to almost 1. This quasiperiodic time dependence from a lower 
stress state to a high stress state is the origin of the nonergodicity for a ^ 1 and R ^> 1. 
That is, small earthquakes accumulate stress locally and characteristic earthquakes release 
the stress globally and quasiperiodically. The periodicit y of large characteristic earthquakes 
for R = 1 and a = 2.5 was also observed in Refs. 20|, |21(. The quasiperiodic behavior in 
Fig. [ED is reminiscent of the stress versus time curves observed in laboratory experiments 
with rocks 



50]. 



For R = 1 and a = 0.5 46j, there are intervals where the stress increases after small 
events. However, many random decreases occur and the time series fluctuates between 0.3 
and 0.8 similar to the behavior for a = 2.5 and R = 1. These random fluctuations are 
consistent with the system being ergodic (see Fig. [9]). In contrast with the behavior of 
the time series for a = 2.5, no quasiperiodic behavior is observed as R is increased (see 



Figs. 11(c) and 11(d)). Instead, f(t) remains in a high stress state, </> ~ 0.99, and the 



system remains ergodic. This behavior is also observed in the long-range CA models 32]. 
Note that for R = 100, the stress returns to a relatively small value only once during 
the observation time, which makes the system exhibit punctuated ergodicity as shown in 
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FIG. 11: Time series of the stress / per block as a function of substrate plate updates (loading 
times) for (a) R = 1, a = 2.5, (b) R = 500, a = 2.5, (c) R = 100, a = 0.5, and (d) R = 500, 
a = 0.5. For R = 1 the stress fluctuates between and 0.4. As R is increased, the time series 
becomes quasiperiodic and the system is no longer ergodic. The sharp decrease in / for R = 100 



and a = 0.5 corresponds to the punctuated ergodicity seen in Fig. 10(b) As R is increased, no 



quasiperiodic behavior is observed, and / fluctuates about a high stress state close to 1 — a. 



Fig. 10(b) 



Increasing a for R = 500 makes the quasiperiodic behavior of f(t) better denned and 
increases the period 46[]. For a = 0.9 (not shown) 46J the system remains in a high stress 
state for some time and exhibits mode- switching behavior similar to the long-range Olami- 
Feder-Christensen model with a long healing time [|5l|], and a cellular automaton model of a 
vertical fault with dynamic weakening of cell strengths 52]. Ben-Zion et al. [53| have noted 



the importance of mode switching in understanding earthquake fault systems. 
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FIG. 12: (Color online) The waiting-time distribution P(t) for (a) R = 1, a = 2.5, (b) R = 500, 
a = 2.5, (c) R = 1, a = 0.5, and (d) R = 500, a = 0.5. The time is scaled by the mean waiting time 
<t>. The waiting-time distributions are fitted to the exponential form f(t) oc e _ */ T if possible. 
The characteristic time r of the exponential decay increases with M for the power law events. 
The waiting-time distribution P(t) is close to an exponential for both the scaling events and the 
characteristic events for R = 1. For R = 500 the waiting-time distribution for the scaling events 
is close to an exponential, but P(t) for events with M > 1 is not. For clarity, each distribution is 
shifted vertically by 0.1. 



VIII. WAITING-TIME DISTRIBUTION 



The nature of the waiting or recurrence-time distribution P(t) for events of a given range 
of sizes is of current interest 54|, |55|, [56] . Statistical data from the Southern California Fault 
Network show there exists correlations, at least between large events 54], |55|, |56[]. These 
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correlations imply that the waiting-time distribution is not exponential. 

We assume that the substrate moves with a constant velocity and hence the time t between 
two events is proportional to the number of substrate updates. Figure [T2] shows our results 
for P(t) for a = 0.5 and 2.5 and R = 1 and 500. For all combinations of R = 1 and a the 
waiting-time distribution for both the scaling events and characteristic events is close to an 
exponential distribution, which implies that there is little correlation between the events. 



The exponential decay of P(t) for R = 1 and a = 2.5 has been reported in Refs. 20j, |21 | 



For R = 500 and a = 2.5 (see Fig. 12(b) ) the waiting time distribution for the characteristic 
events is nonexponential because the characteristic events are quasiperiodic for large R. 
Because there were no characteristic events for R = 500 and a = 0.5 during our observation 
time, we computed P(t) for events in the range 1 < M < 10 and found that there is a 
maximum at t pd 100 and a long exponential tail. Hence, we conclude that the large events 
in long-range models for a = 2.5 and a = 0.5 are correlated. 

IX. SUMMARY 

The Burridge-Knopoff model with long-range stress transfer is more realistic than the 
usual Burridge-Knopoff model with nearest neighbor interactions and is more realistic than 
the long-range CA models because of the presence of inertia and a dynamic friction force. 
Our results show that the generalized Burridge-Knopoff model exhibits much richer scaling 
behavior than the cellular automata models and its behavior depends on the nature of the 
velocity- weakening friction force and the range of the stress transfer \&\ . 

For nearest-neighbor interactions (R = 1) we verified the power law behavior of the 
moment distribution P(M) for M ^ 1.0 and a ^ 1 with a scaling exponent of x ~ 2 for 



a 



> 2 



14 



, 1 15l | . Qualitatively similar results (not shown) 46| were found for a = 3, 4, 5, 
and 10 with £ = 10 and for £ = 5, 7, 10, and 20 with a = 2.5. No scaling was found for the 
moment distribution for a ^ 1, and the distribution of the number of failed blocks during 
an event, P{s), does not scale for all values of a considered. 

For the long-range model (R 1), P(M) and P(s) show similar behavior. For a ^ 1, 
power law behavior with an apparent exponent of x ~ 2 as for R = 1 was found with 
the additional presence of characteristic events. As R is increased, the range of power law 
behavior becomes smaller. For a ^ 1, an exponent of a; ~ 1.5 was obtained for R ^> 1. This 
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value of x is consistent with that found for the long-range CA models and the existence of 



a spinodal critical point 



29. 



30 



31 



32] . The probability of characteristic events for small a 



decreases rapidly with increasing R. 

We found that the Burridge-Knopoff model is ergodic for R = 1 and all values of a 
studied; this behavior is in agreement with the random fluctuations in the time series of the 
stress. For R ^> 1, the system is ergodic for a < 1 because the system fluctuates about 
a high stress state. This behavior is consistent with the behavior of the long-range CA 
models 



30 



31 



321 ] and the observation of the Southern California fault system |48j. The 



system becomes nonergodic for R ^> 1 and a ^ 1 due to the quasiperiodic behavior of the 



stress. For a = 0.9 and R = 500, we found mode-switching between quasiperiodic 



Dehavior 



and fluctuations around a high-stress state, similar to that found in other models 51 



52 



53]. 



The exponential fits to the distribution of waiting times for the scaling events for the values 
of a and R studied implies that there is no correlation between these events, which is 



inconsistent with the statistical data from real fault network systems f 
events such as characteristic events are correlated. 



. However, large 



Our simulation results suggest that there exists two scaling regimes with qualitatively 
different behavior, one of which (a — > and R ^> 1) is consistent with an equilibrium 
spinodal critical point and a mean- field exponent of x = 3/2, similar to the long-range CA 



29 



15 



3G, 31 1 . The nature of the 

n 

, 118|] . The apparent 



models and the near-mean-field picture of spinodal nucleation 
other scaling regime with x ~ 2 for a ^ 1 is not well understood 14, 
dependence of the scaling exponent x on a for 1 < a ^ 2.5 suggests that the interpretation 
of this scaling regime in terms of a dynamical critical point must be viewed with caution 
and that larger system sizes as well as longer run times should be investigated. Because 
real faults are finite and the number of events observed is small, the a-dependence of x seen 
in the long-range Burridge-Knopoff model may accurately reflect the behavior of some real 
faults. The qualitatively different behavior observed for a ^ 1 and a ^ 1 is consistent with 
recent results Q for R = 1. 

For large a and large R, our results resemble those observed in laboratory experiments 
on rocks. As the range R is increased and a is decreased, our results more closely resem- 



ble the long-range CA models. This wide range o 



several models of earthquake faults 



14 



behavior indicates that the physics of 
241 ] can be obtained from the generalized 



Burridge-Knopoff model with the appropriate choice of R and a. Our results can be in- 
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terpreted as suggesting that real earthquake faults and laboratory rocks can have different 
statistical distributions of events and different physical characteristics due to the details of 
the friction force as well as range of the stress transfer. This dependence means that we must 
develop ways of determining the friction force with considerably more accuracy if we want 
to understand the relation between physical processes and observed earthquake phenomena. 
Our waiting-time results also shows that other physical features might need to be added to 
the general^ Bnrridge-Knopoff mo del considered in this paper su,h „ ry effects, 



heterogeneities, and a more realistic rate-dependent friction force [21. l59l. l60l|. 
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